The random K-satisfiability problem: from an analytic solution to an efficient 

algorithm 



Marc Mezard^ and Riccardo Zecchina^ 
^ Laboratoire de Physique Theorique et Modeles Statistiques, 
- - - CNRS and Universite Paris Sud, 

Bat. 100, 91405 Orsay CEDEX, France 
»' , ^ The Abdus Salam International Centre for Theoretical Physics, 

^ ' Statistical Mechanics and Interdisciplinary Applications Group, 

Str. Costiera 11, 3^100 Trieste, Italy 
^ ; (Dated: February 1, 2008) 

►-5 . We study the problem of satisfiability of randomly chosen clauses, each with K Boolean variables. 

' Using the cavity method at zero temperature, we find the phase diagram for the K=3 case. We 

I show the existence of an intermediate phase in the satisfiable region, where the proliferation of 

metastable states is at the origin of the slowdown of search algorithms. The fundamental order 
parameter introduced in the cavity method, which consists of surveys of local magnetic fields in the 
various possible states of the system, can be computed for one given sample. These surveys can be 
^ , used to invent new types of algorithms for solving hard combinatorial optimizations problems. One 

H ■ such algorithm is shown here for the 3-sat problem, with very good performances. 
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I. INTRODUCTION 



. The K-satisfiability (Ksat) problem deals with an ensemble of N Boolean variables, submitted to M 
^ constraints. Each constraint is in the form of an 'OR' function of K variables in the ensemble (or their 
^ negations), and the problem is to know whether there exists one configuration of the variables (among 
I I - the 2^ possible ones) which satisfies all constraints. The Ksat problem for iiT > 3 is a central problem 
\ in combinatorial optimization: it was the first problem to be shown NP-complete ^, and an efficient 
. algorithm for solving Ksat in its worst case instances would immediately lead to other algorithms for 
^ ' solving efficiently thousands of different hard combinatorial problems. 

At the core of the statistical physics of disordered systems is the spin glass problem (SG), which 
also deals with Boolean variables ('spins'), interacting with random exchange couplings |^. Each pair 
, of interacting spins can be seen as a constraint, and finding the state of minimal energy in a spin glass 
' amounts to minimizing the number of violated constraints. Although the precise form of the constraints in 
C^l SG and Ksat differ, there exist deep similarities ||]; in both cases the difficulty comes from the existence 
of 'frustration'!^, which forbids to ffiid the global optimal state by a purely local optimization procedure. 
Links between combinatorial optimization and statistical physics have been known for longQ. Two main 
categories of questions can be addressed. One type is algorithmic, for instance finding an algorithm which 
decides whether an instance is SAT or UNSAT. Another is more theoretical, and deals with large random 
^ instances, for which one wants to predict the typical behaviour. Examples of use of statistical physics 
(— ( in each category are the simulated annealing algorithm Q and the solution of the random assignment 
Q \ problem [0, ^ , or the direct mapping of certain graph partitioning problems to spin glasses |^ . Here we 
O ■ address the two types of questions in the 3sat problem. 

^ ' The study of random Ksat problems, where the clauses are chosen randomly, is also interesting from the 
viewpoint of optimization. In practice, algorithms which are used to solve real-world NP-complete prob- 
lems display a huge variability of running times, ranging from linear to exponential, when the parameters 
(e.g. the number of clauses) are changed. A theory for the typical-case behaviour of algorithms, on classes 
of random instances chosen from a giv en p robability distribution, is therefore the natural complement to 
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the worst-case analysis |1( 

In random 3sat, numerical simulations have shown the existence of a phase transition when one varies 
the ratio a = M/N of the number of clauses to the number of variables. For a < ac the generic problem 
is satisfiable (SAT), for a > ac the generic problem is not satisfiable (UNSAT) Using the cavity 
method, first developed in spin glass theory, we shall show the existence of this threshold and compute 
ttc — 4.267. We also find an intermediate region, the HARD-SAT-Phase ad = 3.921 < a < ac, where the 
generic problem is still SAT but the proliferation of metastable states makes it difficult for algorithms to 
find a solution. This proliferation is similar to the effect found in the theories of structural glasses using 
spin glass models with multispin interactions RTj , where it is known to lead to a dramatic slowdown of the 
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relaxation. In this sense the difficulty to solve the 3sat problem in the intermediate region ad < a < ac 
is similar to the difficulty in equilibrating structural glasses. 

This theoretical analysis is done using the cavity method, at a level equivalent to what is called one 
step replica symmetry breaking in the replica language. This means that it assumes the existence of many 
states, but cannot handle a nontrivial correlation pattern between them. There are some arguments which 
point towards the correcteness of this solution, although an exact proof looks somewhat remote at present. 

In this cavity method with many states, the order parameter consists in the surveys of local magnetic 
fields acting on each spin. While for the theoretical analysis one averages over the random graph structure 
of the problem, it turns out that this order parameter can also be computed for one given sample, using 
a reasonably simple message passing procedure which takes into account the multiplicity of states. This 
procedure provides a generalization of the belief propagation used in statistical inference ; it is shown 
here to converge in some difficult situation with many equilibrium states where ordinary belief propagation 
does not converge. The resulting surveys provide an interesting description of one given sample, where 
the various variables are found to play very different roles. This single sample analysis is very useful in 
order to find new algorithms for solving hard optimization problems. Here we show one such algorithm 
for the random 3sat problem, where the surveys are used to identify one spin and fix it. The problem is 
thus reduced and the surveys are computed again on the new system. This decimation procedure is shown 
to have very good performances, comparable to, or better than, the state of the art in this problem. 

The paper presents a number of concepts and techniques, both analytical and numerical, which can be 
applied to a rather large class of combinatorial optimization problems. We have presented these concepts 
and techniques in a broad framework, in order to allow for future use on different problems. The concrete 
implementation is then done on the random 3sat problem. Some of the results discussed in this paper 
have been recently announced in | [T^ . 

The paper is organised as follows: in sect. || we present the generic structure of the optimization 
problems in which we are interested. These can be represented as bipartite graphs called factor graphs. 
Sect. Ill recalls a general message passing procedure which can be used to study optimization or inference 
problems defined on these factor graphs. The basic ingredients of this procedure are messages which 



we call cavity-biases which play a crucial role in the whole paper. Sect. |V| defines the set of random 
graphs which we study, which are random hypergraphs with a fixed connectivity. Sect. ^ provides some 
background on the decomposition of the configuration space into states. It consists of a short introduction 
for non-physicists, a specific definition of zero temperature states for the random 3sat problem, and 



the definition of the complexity which is a crucial concept for a system with many states. Sect. VI 



provides an introduction to the zero temperature cavity method, presented in the general setting of 
combinatorial problem on random factor graphs. This section summarizes, and puts in a more general 
context, some recent work which has developed the cavity method for finite connectivity problems, first 



at finite temperature |l9[| , then at zero temperature 1 20 . It provides the whole formalism for t he a nalytic 



study of the phase diagram. This formalism is applied to the random 3sat problem in sect. VII , where 
all the results o n the phase diagram are derived. We explain the survey propagation algorithm on a given 
sample in sect. VIII , and the decimation algorithm for solving large random 3sat problems is presented 
in IX. Sect. ^ contains some concluding remarks. A ppen dix A contains some technical details relative 
to the computation of the phase diagram done in sect. VII. Appendix B explains the computation of the 
free energy for one given sample. 



II. FACTOR GRAPH REPRESENTATION 



The models we are interested in involve Boolean variables which interact in groups, the energy being 
the sum of energies of all groups. We shall adopt the factor graph representation [2l] familiar in computer 
science, but we shall keep to the representation of Boolean variables as Ising spins, more familiar to 
statistical physicists. 

We consider a set of N Ising spins, Ci g {±1}, and we suppose that we have M groups of interacting 
variables, which are called function nodes. Each function node a involves a set of spins. We denote by 
Va the set of all these spins. The interaction is an arbitrary function of the spins in Va, which depends on 
the problem one considers, and can also involve hidden variables. 

The total energy of a configuration ai, ■■■,(Jn is 



M 

a=l 



(1) 
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FIG. 1: An example of a factor graph with 5 variable nodes 4 = 1, .., 5 and 3 function nodes a = 1,2, 3. In this 
case, each function node has connectivity 3, as in the 3sat problem. The connectivities of the 5 variable nodes are 
respectively 2, 2, 1, 1, 3. 

and the goal in combinatorial optimization is to find a configuration of spins which minimizes E. A 
generalization of this problem, natural from the point of view of physics, and which connects with problems 
in statistical inference, consists in introducing an additional parameter /3, an 'inverse temperature' in the 
physics language, and in studying the Boltzmann probability distribution: 

P{ai,...,aN) = |exp(-/3S) , (2) 

where Z is a normalisation constant. As usual in physics, we shall denote by < O > the expectation of 
an observable O (which can be any function of the ai) with respect to this measure. In the large (3 ('low 
temperature') limit this measure concentrates onto the lowest energy configurations. At finite (3 one may 
be interested in computing for instance the expectation value of one spin variable tri with the Boltzmann 
probability. Because we work with binary spins, this average determines the full marginal probability law 
of CTi: this is precisely the quantity that one typically seeks in many inference problems, like e.g. decoding 
procedures for error correcting codes. 

The general problem can be represented by a graph consisting of two types of vertices, 'variable nodes' 
associated with each spin, and 'function nodes'. A function node a is connected by edges to all the variable 
nodes involved in Sa- Therefore each variable node has connections towards all the function nodes in which 
it appears, and the graph is bipartite (see fig.|l|). Each spin ai is connected to rii function nodes, we denote 
this set of function nodes by Vi. We call n.i the connectivity of spin i, Ua the connectivity of function 
node a. Throughout this paper, the variable nodes indices are taken in k, while the function nodes 
indices are taken in a,b,c, ... 

Let us give here a few standard examples. 
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Spin glasses: all the interactions involve two spins, so all Ua are equal to 2; the energy of an interaction 
node a involving spins ct.; and aj is given by 

Ea = -~Jij(Ji(Jj , (3) 

where the number Jij is called the coupling constant. Generalized spin glasses with p-spin interactions, 

Ea ^ -Jii...i^cri^...cri^ , (4) 

have also been studied a lot in statistical physics as models of structural glasses. They are the closest 
physical analogues of the satisfiability problems which we study here p2| . 

K-sat: All interactions involve K spins, and the energy of an interaction node a involving spins 
CT/i, ...,cri^ is given by: 

E. = 2\[^l±^. (5) 

r=l 

It depends on a set of K coupling constants Ja = {J^, Ja) which take values ±1. This interaction 
node has a simple interpretation as a clause: the energy Ea is zero as soon as at least one of the spins 
is opposite to the corresponding coupling J^. If all spins are equal to their couplings, the energy is 
equal to 2. The more conventional description of Ksat uses Boolean variables: let us introduce Xi which 
is TRUE if and only if ai — 1. The energy Ea depends on the OR of the K variables yi^ , where yi^ 
is the original Xi^ when = —1 and is its negation when = +1. The energy vanishes if j/i^ V ... V y^^^ 
is true (the clause is then said to be satisfied), otherwise it is equal to 2 and the clause is unsatisfied. This 
arbitrary factor of 2 is introduced for future convenience. 

One can also consider graphs involving mixtures of function nodes of different types, e.g. mixtures of 
K — 2 clauses and K = 3 clauses |2^. These are some examples of constraints satisfaction problems, 
but of course there exist many other instances of problems, much studied in computer science, which can 
be represented by such factor graphs. 

In general, an instance of the problem (also called a sample in physics language) is given by a graph and 
the set of couplings needed to define each function node. In physics (e.g. spin glasses) one is interested 
in the configurations of low energy. In the SAT problem, one wants to know whether there exists a 
configuration of zero energy (in which case the instance is called SAT), or not (in which case the instance 
is UNSAT). 



III. THE SUM-PRODUCT ALGORITHM 

A popular method for studying the inference, i.e. the probability measure (H), is a message passing 
procedure called the "Sum-Product" algorithm [pT| When used at (3 ^ oo, the corresponding "Min- 
Sum" algorithm can also be used to get some information on the lowest energy configurations. This 
procedure is exact and fast on tree- like graphs. In our case the Sum-Product algorithm amounts to 
sending some messages along the edges of the graph. We call cavity-field, and denote by hi-ta, the 
message passed from a variable node i to a function node a. We call cavity-bias, and denote by Ua^i, the 
message passed from a function node a to a variable node i. 

The cavity-field hi^a is given by the sum of cavity-biases converging to i from all function nodes b 
distinct from a: 

hi^a = ^ Ub^i (6) 
beV{i)/a 

The operation performed by a function node to compute the cavity-biases which it will send to its 
neighboring variable nodes is a partial summation: it computes the marginal probability law for that 
variable to which it sends the message. More precisely, let us consider a function node a of connectivity K 
and let us suppose, for notational simplicity, that it is connected to variables cti, uk- The cavity-bias 
Ua^i sent from the function node a to the variable node i = 1 is a function of the cavity-fields hj^a sent 
from all other variables nodes j E {2,...,K} towards node a. One considers the function of aidefined 
by J2a2,...,crK exp(-/3£^Q(cri,CT2, ...,(7^) + P {h2->aO-2 + ■■■ + hx^aCTK))- As (Ti = ±1, this fuuctiou can be 
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written for instance as the exponential of a linear form in cri. The cavity-bias Ua^i sent from a to 1 is 
defined from: 

exp{-/3Ea{ai,a2,--,(TK) + P {h2^a(J2 + ■■■ + hK^aCTK)) = exp{f3 {Wa~,l + (JiUa^l)) ■ (7) 

cr2,...,crK 

Beside the cavity-bias Ua^i, this equation also defines a 'free-energy shift' Wa—>i which is not used in 
the Sum-Product algorithm but will become very important in our generalisation later on. In physics 
words, h.i_^a is the magnetic field on spin number i whenever the interaction a is turned off, and Ua^i is 
the contribution to the magnetic field on spin number i from the interaction a. Equation ^) indicates 
that the probability law of spin ai due to the interactions a = 2, ..K is a product of independent laws 
due to each interaction, while the marginalization operation (|^) is a partial summation, hence the name 
sum-product. The algorithm is easily generalized to variables which are more complicated than Boolean. 

The iteration of the above message passing algorithm, starting from a generic random initial condition, 
is known to converge whenever the underlying factor graph is a tree. Actually it converges in one sweep 
if one first computes the messages from the leaves of the tree. The resulting set of messages can be used 
to compute the probability distribution of one spin (or more generally of some subset of spins). One just 
needs to compute the local field Hi on spin ai: 

H, = V Ua^^ , (8) 



aeV{i) 



and the probability distribution of ai is: 



_ exp(/?g.a,) 

^^^^^ " 2cosh(/3ifO ■ ^ ' 

One way to prove this result, using the physics language, is to show that the message passing algorithm 
minimizes the Bethe free energy of the spin system [Q. As the Bethe free energy is exact for tree like 
graphs, this provides the proof. If one is interested in the optimization problem (/3 — > oo), one can 
show that the configuration ai = SignHi is the lowest energy configuration if there is a unique such 
configuration. If there are several lowest energy configurations, taking the /? — > oo limit of < ai > with 
the measure gives the average of the spin ai over all these configurations (but one needs to make this 
detour through the finite (3 problem in order to get this result). 

In this work we shall mainly be interested in the optimization problem. Working directly at /3 = oo 
then simplifies the algorithm. The cavity-fields are generated as before by sums of cavity-biases. For 
computing a cavity-bias, one performs a partial minimization, and the formula (0) simplifies to 

min [Ea {ai, a2,. ak) - {h2^aa-2 + ■■■ + hk^aO-k)] = - {wa^i + <yiUa~,i) ■ (10) 

This equation defines the output messages w and u as functions of the input messages hj^a- In general, 
we shall write 

Wa-,1 = Wj^{h2^a, --.hK^a) ] Ua^l = Uj^{h2^a, flK^a) , (H) 

which defines the functions wj^ and iij^ (the label Ja is here to explicitly remind that a given function 
node energy Ea will in general depend on some set of couplings - see (0,^)- which we denote collectively 

as Ja). 



IV. RANDOM GRAPHS AND THERMODYNAMIC LIMIT 



A. Definition of random graphs 



In the rest of this paper we shall consider the random K-sat problem, which is defined on some ensemble 
of random graphs which we now describe. To lighten the notation, we concentrate on the K = 3 cases 
for which all the function nodes have connectivity ria ~ 3, and we generate the random graphs as follows: 
For each triplet i < j < fc of variable nodes, a function node connecting them is present with a probability 
6a /N'^, and it is absent with probability I — 6a /N^. The average number of function nodes is then 
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M = aN . The graph model used is analogous to the G(N,p) model of random graph theory (see, e.g. 
0), withp = 6a/7V2. 

For the problem which we consider, the energy Ea associated with a function node a also depends 
on some coupling constants Ja (see for instance (g|)), which may be drawn randomly and independently 
for each function node from some probability distribution. For instance in random 3sat, each number 
Jq, Jq, takes values ±1 with probability 1/2. In general we shall denote by SjO the average of any 
quantity O over all the choices of the random graphs (with fixed N and M), and over the choices of 
couplings. In such a probabilistic setting, one is interested for instance in computing the average 'ground 
state energy': For each sample, an optimal configuration (one with minimal energy), is called the ground 
state, its energy is Eq, and one would like to compute £jEq. 

B. Thermodynamic limit 

We shall be interested in the Hhermodynamic limif where M and iV go to oo, keeping the ratio a = 
M/N fixed. The connectivities of variable nodes become independent identically distributed (iid) random 
variables with a Poisson distribution fsaik) of mean 3a, since the probability of having k edges connected 
to a variable node is: 

lim f ~ ~ ^ {6a/N'ni ~ ea/TV^)^^-!)^^-^)/^-'^ = ^ exp(-3a) ^ /.^(k) (12) 

N—*oo y rC J fc! 

The structure of the random graphs generated by this process for large TV is interesting. Locally such a 
graph is tree-like: the typical size of a loop in the graph scales like log(7V) for large A''. On the other hand 
loops are definitely present, and they can induce 'frustration' in the sense of having competing constraints 
(it has been argued that similar random graphs with a local tree-like structure provide a natural setting 
for discussing the 'Bethe approximation' of frustrated systems (2^). The structure of the graph has one 
important consequence: consider one given function node, connected to three variable nodes (spins). If 
one deletes this function node, the typical distance between any two of these three spins (measured as the 
length of the shortest path on the graph which connects them) is of order logiV, and thus diverges in the 
thermodynamic limit: the spins are far apart. This property will be crucial in understanding the type of 
correlations existing between the spins, and in solving the model. Notice that the limit where M, iV — > oo 
is also the one that is interesting from a computational complexity point of view. 

For problems defined on random graphs with given N, M, the ground state energy fluctuates from sample 
to sample. It is often true, but it may be difficult to show, that the distribution of the ground state 'energy 
density' Eq/N becomes more and more peaked when N increases, so that, in the thermodynamic limit, 
almost all samples have the same energy density, which can be computed as 

£0 = lim £jEo/N (13) 

For random Ksat it can be proved that the above condition holds . 

One of our aims is to compute this limiting value eo for a fixed value of a = M/N. For the Ksat problem 
it turns out that this value is equal to zero below a certain threshold ac, and becomes > for a > ac- In 
statistical physics it is very difficult to go beyond the estimate of the energy density: if one can compute 
£0, one knows that Eq e^N is the leading behaviour of the energy for large N, but in general one cannot 
control the subleading part, and in principle it could be possible for instance that also in the small a 
phase of Ksat where — 0, some finite contribution to Eq (finite when N oo) could make the problem 
typically UNSAT. Numerical simulations tend to show that this is not the case. Knowing eo will then 
allow to get the phase diagram of the problem. But one is also interested in other properties of the generic 
samples in the thermodynamic limit, like the decomposition of the space of accessible configurations at a 
given energy, to which we now turn. 

V. STATES AND CLUSTERING PROPERTY 
A. A simple example of pure states: the ferromagnet 

One of the main aims of statistical physics is to understand the building up of correlations between 
distant variables, when the basic interactions between them are short range. This is precisely the type of 
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question that we need to address here: variables interact locally (the only direct interactions involve spins 
connected to the same function node). But we also need to control the correlations established between 
two spins belonging to the same function node, due to their indirect coupling through other nodes. As we 
saw, the geometry is such that this indirect interaction builds up through very long (O (log TV)) paths. 

Usually, in the statistical physics of systems with short range interactions, the correlation between 
distant variables displays a relatively simple variety of behaviours. The simplest one is when there is 
only one pure state in the system (typically a 'paramagnetic phase'): then there exists a finite correlation 
length and the connected correlation function between two distant spins , aj decays exponentially with 
the distance dij at large distances: 

I < UiCJj > - <(Ji X o-j > I ~ C exp{-dij/^) (14) 

where ^ is the correlation length (C can be a constant, or involve power law corrections in the distance). 
This is called the clustering property. On the other hand, some systems can also have phase transitions, 
and display a low temperatue phase with several pure state. 

The archetypical case which we briefly describe here as a pedagogical example is the ferromagnetic 
p — 2 spin system with energy given by (^ with Jij = 1: at low temperature the spins polarize in one 
of two 'pure states', related to each other by the global symmetry changing all ai to —ai. Let us call a 
configuration one assignment of the A'' spins, di, aN- The pure states are probability measures on the 
configuration space obtained using a slightly modified Boltzmann measure where one adds an external 
'symmetry breaking' magnetic field (in the present language one adds a function node of connectivity one 
connected to each variable node i , with energy —Bui). One computes lim^^gi liniAf^oo < >i which 
defines the expectation value < Ui >± of spin i in each of the two states + and — . It is a well known fact 
that the connected correlation functions within each state have the clustering property: 

I < ai(jj >± - < cr,; >±< cTj >± I ~ C cxp(-dy7^) (15) 

This means that when the spins collectively polarize in the + state, the correlations between distant spins 
vanishes. It is not true for the full Boltzmann measure: if one does not add the small symmetry breaking 
field B, but keeps to the Boltzmann measure (^, one gets for any observable < O >— l/2(< O >+ + < 
O >-) (the fact that the two states enter with equal weight 1/2 is a consequence of the global symmetry 
of the original problem), and one easily shows that the corresponding correlations do not vanish at large 
distances. 

The lesson we learn from statistical physics is that correlations decay at large distance within each pure 
state. In problems more complicated than the ferromagnet it may be difficult to identify the various pure 
states, especially when we do not have at hand a simple breaking of a symmetry. A large part of the work 
on spin glasses has been devoted to this problem and we shall not try to reproduce it here (see for a 
review), nor to give a general definition of states at finite temperature in our problems. 

B. States at zero temperature 

Instead we shall focus on the zero temperature limit (/3 — s- oo), where the situation is simpler. A state is 
defined in the thermodynamic limit as a cluster of configurations, all of equal energy, related to each other 
by single spin flip moves, and which are 'locally stable', in the sense that the energy cannot be decreased 
by any flip of a finite number of spins. In the random 3sat problem one can use an even simpler definition 
which allows to generalize the definition of states also to finite N problems: the condition of local stability 
can be substituted by a condition of stability with respect to any sequence of one spin fiips. The reason for 
this simplification, specific to the Ksat problem, is that in this case the stability with respect to sequences 
of single spin fiips insure stability with respect to collective fiips of finite sets of spins. 

C. Many states: definition of the complexity 

Experience with disordered and frustrated systems like glasses shows that there can exist many states, 
and the number of states typically grows exponentially with the number of variables. The number of 
states Af{E) with energy E is written as: 



Af{E) = exp{Ni:{a,e)) 



(16) 
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where the quantity T,{a,e) is cahed the complexity. It is a function of a = M/N and e = E/N, and 
the form (|l^) is derived from the basic assumption that \ogM{E) is extensive. In general, whenever a 
problem has a nonzero complexity, one may expect that simple local algorithms will have great difficulty 
in finding the ground state, simply because the states proliferate (for large TV) and the algorithm will 
easily get trapped into one state with energy above that of the ground state. We shall see in the next 
sections how the cavity method can handle such a situation. 

VI. A PRIMER ON THE CAVITY METHOD AT ZERO TEMPERATURE 

The cavity method was originally introduced in [p6| to study spin glasses, but it gives a general frame- 
work for computing statistical properties of various frustrated systems, and is ideally adapted to systems 
with a locally tree-like structure. It is always in principle equivalent to the replica method, which is 
a more compact and very appealing formalism, however it possesses two advantages. On one hand, it 
proceeds through a standard probabilistic analysis, and makes explicit all the hypotheses involved in it. 
Roughly speaking, the cavity method assumes some properties about the correlations between variables 
in a system with N spins, and shows that these are self-consistently reproduced for a system with iV -f 1 
spin system. The problem in turning it into a rigorous methods is that these hypotheses only hold in the 
large N limit, not for N small. If one is able to have a good control of the correlations as a function of 
N , then the cavity method becomes also a choice method for rigorous probabilistic studies of frustrated 
systems On the other hand, in the cavity method, one considers explicitely the site dependence of 
the order parameter, and the averaging over 'disorder' is performed at the end (this is in contrast with 
the replica approach where the disorder average is made from the very beginning). As we shall see here, 
this aspect allows to define some algorithm, inspired from the cavity method, which computes the order 
parameter on each site for one given sample. 

In what follows we shall be interested in the zero temperature version of the cavity method. As discussed 
in details in ||2^, the formalism simplifies a lot in this limit. Here we shall mainly outline for completeness 
the basic aspects of the method, applied to the 3-sat model where all function nodes involve exactly three 
spins (the generalization to more general problems is totally straightforward but would make the notation 
more cumbersome). We refer the interested reader to for more details. We shall first present the 
method in its simple 'replica symmetric' (RS) version where it assumes the presence of a single state, and 
we shall then turn to the more involved case in which many states exist but are uncorrelated, a situation 
called 'one step replica symmetry breaking' (IRSB) in the replica jargon. 

A. The cavity method with one single state ('RS' case) 

1. Adding one spin 

Consider a N spin system ai,..,aN and its interaction graph, and add to it a new spin cto. Then 
generate the new function nodes involving this new spin as follows: for each pair 1 < i < j < N , the 
function node (0,i,j) is present with probability 6a/N^. Therefore we have added k new function nodes 
which we label by a = 1, fc, where A: is a random variable with probability distribution fsaik). Let us 
consider all the new function nodes which involves, besides (Tq, 2k other spins which we call and 
(see fig]^). Generically, on the original graph (i.e. before adding ctq), these spins are far apart from each 
other. If there exists only one state, the clustering property implies that the correlations between these 
spins, before adding ctoj vanish. Using the fact that and are binary variables, this decorrelation 
implies that the minimal energy of the original graph, for fixed values of the 2k spins and cr^, can be 
written as 

k 

E{{alal})^A-Y,{gaal + h,al) (17) 

a=l 

where the 2k 'local fields' Qa and ha are nothing but cavity-fields passed from each spin to the function 
node a, and ^ is a constant (independent of the local fields). 

Looking at the function node a in the full graph including spin cto, we need to minimize the function 
Eaio'OjO'ai'^a) ~ (Sa'^a + ^aO'a) with rcspcct to o'Q,cr^. This is precisely what one does in the message 
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passing procedure of sect. [11, and one can thus use (jT^) to get the minimal energy of the new graph with 
iV + 1 spins, for a given value of ao: 



E{aQ) = ^ - XI ^Ja(ffa, ha) - O-O X UlAga, K) 



(18) 



a=l 



a=l 



where ~ A — X]a=i (l^^al + I'^ol) the minimal energy of the N spin system. 

Equation ( p^ ) shows that the 'cavity field' on the new spin ctq (the coefficient of — ao) can be written 
as: 



ho ^^uj^{ga,ha) . 



(19) 



As we shall see, its is often useful to decompose this cavity-field as a sum of cavity-biases: 

k 



ho = '^Ua ; Ua = UlSgai ha) 



(20) 



2. Self- consistency equation: the order parameter 

Whenever one adds a new spin ctq, one picks up a value of k, and a set of 2fc fields ga,ha, which are 
iid variables taken from a probability distribution V{h). The cavity method assumes the existence of a 
thermodynamic limit iV — s- cx) where the energy density E/N and the distribution of local fields V{h) have 
well defined limits. This means that the distribution of Hq is the same as that of the 2k fields. Calling 
Q{u) the probability distribution of the u variables, this stability condition of the iteration(|2^) implies 
that: 

Qiu) = j dgdhV{g)V{h)£j5{u-u3{g,h)) 

- E/3«(fc) / dui...duuQ{ui)...Q{uk)5\h-Y,Ua\ (21) 

fe=0 \ a=l / 

where £j means an expectation value with respect to all the couplings J . 
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3. Computing the energy 



One can easily compute the average s hift in the ground state energy when adding one new spin. Looking 



is A-J2a=ii\9a \ + \ha\), while that of the iV + 1 spin system is ^- X]a=i ^j(3a, ha)-] J2l=i 
Therefore the shift in energy when adding the new spin ctq is 

k k 

Ae[°^ ^J2'^-'^J^(9a,K) + \9a\ + \ha\)-\Y.ujA9a,ha)\ . (22) 

a— 1 a— 1 

Equation ( pl| ) gives an integral equation for the 'order parameter' which is the probability distribution 
Vih) (or alternatively Q{u)). Let us now suppose that this equation has been solved (we shall see below 
how this can be done on each specific example), and show how the energy density can be deduced from 
this order parameter. We must compute the average of the energy shift for adding one spin (averaged 
over the choice of k and of the corresponding cavity- fields): 

AEi = f2h4k)£j I \{[dgaV{ga)dhaV{ha)] 

k=0 •' a=l 

/ k k \ 

[-*Ja(.9a, ha) + \ga\ + \ha\] - \ ^JaCffa, ha)\ (23) 



at the addition process defined in sect. VI A 1, we see that the energy of the original graph with N spins 



k , , \ I v^fc 



One might believe that, as the energy grows linearly in N at large A^, this average energy shift would be 
equal to the energy density; however there is a correction term due the change in the number of function 
nodes per variable in the iteration N ^ N -\-\. Indeed in the N ~\-\ spin system we are generating function 
nodes with probability Qa/N"^ in a system with 1 vertices and therefore we are slightly over-generating 
function nodes. We need to cancel a fraction 1 — N"^ / {N + 1)^ ~ 2/A^ of them at random: the probability 
of deleting k' function nodes is 

^ {2/N)''{l ~ 2/iV)"^-'^' ^ (2a)'=' exp(-2a)/fc'! = ha{k') (24) 

and the average number of deleted function nodes is 2a. Each deleted function node contributes to the 
average energy change with a correction term: 



AE2 ^ £j j dhidh2dh3Vihi)Vih2)Pihs) 



min [£'(cri,cr2,o'3) - /iiCTi - /i20-2 - /iso-s] + + |^2| + l^sl (25) 

fi ,""2 / 

The ground state energy density is finally given by: 

eo = AEi - 2aAE2 (26) 



B. The cavity method with many states ('one step rsb') 



1. Iteration within one state 



Let us now see how the cavity method can be used to handle a situation in which there exist many 
states. As far as the clustering condition holds within each state, the iterative method can still be applied 
to each state. The problem is that the iteration induces some crossings of the energies of the states, 
and one needs to take this effect into account properly. We proceed as in the previous section by adding 
the new spin cto connected to the 2k spins {o';J,o'^}- In each state a, one can reproduce the previous 
arguments: due to the vanishing of correlations, the energy of the state a, for fixed values of the 2k spins 
{cr^} and {cr^}, can be written as 

k 

E'-iWlal}) = - E i9>l + hyl) . (27) 

a=l 



11 



We now have, for each state a, 2k 'local fields'. 

Withm each state a, the optimization procedure on the 2k spins cr^, cr^ proceeds as before. The minimal 
energy of the new graph with + 1 spins, for a given value of (Jq, is: 

fe k 

E"{ao) = A" ~-Y.'^j^(9^,K) - aoY.ujA9:,K) ■ (28) 

a—1 a—1 

This shows that the local field on the new spin (Tq in state " can be written as 

k 

h^^Y.^jA9a,K) , (29) 

a=l 

and the shift in energy of this state is (see (p2|)): 

k k 

AE^ = J2 i-^jAg:, K) + K\ + \K\) - I E K)\ ■ m 

a—1 a—1 



2. Hypotheses on the states; u-surveys 

We suppose the existence of many states, with a complexity function S(q;, E/N) defined as in (|l6|) which 
is an increasing convex function. Let us consider all the states a with a given energy density E/N — e. We 
suppose that all the local fields h"^^ on a given edge j ^ a are iid, taken from a probability distribution 
Pj^ai.^) called a h-survey. This probability distribution fluctuates from one edge to the next, so that 
the full order parameter, obtained by averaging over edges, is the functional probability distributions of 
these h-surveys. The same hypotheses hold for the distribution of the cavity-biases: all the Ua^o on 
a given link are iid taken from a probability distribution Qq^q(u) called a u-survey. Notice that the 
previous 'RS' solution corresponds to having deterministic messages on each edge Pj^^W — 5{h — hj-,a) 
and Qq^o(u) — 5{u — Ua-to). In the many state hypothesis, we have to generalize the messages, and the 
important quantities are the probability distributions (over the many states having a fixed energy density) 
of the cavity-biases going through a given link. 



3. Iteration: level crossings and reweighting 



One iteration step of the cavity procedure leads to an equation relating the probability distributions, 
before any averaging over the graph. In our iteration procedure, the h-survey on the new site, Pq (/i), is 
related to the h-surveys on the other 2k spins o\,a^ computed in the absence of ctq- Let us denote by 
Paiga) the h-survey incoming onto and P^'{ha) the h-survey incoming onto a^. The h-survey P^ih) 
is given by: 



k 

PSih) = C H [P!iga)dga P!'iK)dK] S 

a=l 



exp 



a=l 
k 



y^iw3Aga,ha) \ga\ \h-a 

a—1 a—1 



(31) 



where C is a normalisation constant insuring that Pq (/i) is a normalised probability distribution. This 
equation provides the generalization to the RSB case of the simple iteration ( [20| ) of the previous section. 
Two complications have appeared: the simple messages (cavity-fields and cavity-biases) have become 
h-surveys, i.e. probability distributions of simple messages, and a new term has appeared which is the 
exponential reweighting term. In this term, the parameter y is a number equal to the derivative of the 
complexity with respect to the energy. 
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Let us now explain the origin of this new and crucial term. For a given state a, we add one new spin 
to the system and want to compute the new h-survey. In this process, we have seen in (|3^) that there 
is an energy shift Ai?" which depends on the state, and is correlated to the value of the cavity-field 
h — X)a=i ^j(5a J )• Let US Call SQ{h,AE) the joint probability (when looking at all states) of the 
cavity-bias and the free energy shift, for this function node. 

When we compute Po{u) at a fixed energy density e = E/N, we get a contribution from all states with 
energies before iteration equal to E — AE. Therefore: 

P^{h) = C J d{6E)So{h,AE)exp (^iVS ~ ^' J '^i^^)^o{h, AE) exp{-yAE) . (33) 

The reweighting term in exp(— yA_B) is due to the level crossing, and the fact that the complexity E(a, e) 
is not constant, but increasing: states with a negative value of the energy shift are thus favoured. 

It is useful, as before, to decompose the iteration procedure (|3l| ) into two steps and introduce the u- 
surveys. On any function node a, we merge two h-surveys Pa(ga) and P^' [ha] in order to build a u-survey: 



Ql{u) = J dgdhP!{g)P:'{h)S {u ~ u,^{g, h)) exp {y [z&jjg, h) - \g\ - M) (34) 
Then we can combine all the u-surveys incoming onto the new spin in order to build its h-survey: 

P^{h)^ J dui...dukQt{ui)...Ql{uk) exp(^y\Y,Ua\^ s(^h-Y,u}j . (35) 

Note that there is some degree of arbitrariness in the way one distributes the reweighting between the two 
iteration steps: different choices amount to different definitions of u-surveys. The above one is the most 
natural one, and this is what we shall adopt from now on. 



4- Order parameter and self-consistency: population dynamics 

Equations (^,^5|) are the main result giving the way to compute the messages sent along to a new 
added site. The assumed existence of a thermodynamic limit allows in principle to write a self consistency 
of the iteration in a way similar to (|2l|). In the present case, this is an equation for the functional V{P{h)) 
giving the probability, when one picks up an edge i — > a at random, to observe on this edge a h-survey 
Pi^a{h) equal to P{h). Alternatively, one can use the functional Q{Q{u)) giving the probability, when 
one picks up an edge a ^ j at random, to observe on this edge a u-survcy Qa-tj{u) equal to Q{u). In the 
following we shall rather work with the u-surveys which turn out to have a simpler structure in practice, 
but obviously a fully equivalent description can be obtained working with h-surveys. 

These functional equations are the generalisation to the RSB case of the equations ( ^l| ) for the RS 
case. One could write them explicitely, but they are not particularly illuminating, and we prefer to work 
directly with the iteration equations ( ^4p^ ). These define a stochastic process; at each iteration, one 
performs the following operations: 

• Pick up at random a number of neighbours k, with the probability fsaik); 

• Pick up at random k u-surveys Qi{ui), ...,Qk{uk) from the distribution Q{Q{u)); 

• Compute a h-survey Pi{g) as the reweighted convolution 

Piig)^Ci y (iwi...dufc(3i(Mi)...(5fc(ufc)exp ^Ua|^ ^ (^S ""^^u^ . (36) 

• Pick up at random a number of neighbours k' , with the probability /3q(/c'); 

• Pick up at random k' u-surveys Qk+iiui), ...,Qk+k'{uk') from the distribution Q{Q{u))] 
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• Compute a h-survey P2{h) as the convolution 

P2(/i) = C'2y"dMi...dufcQfc+i(wi)...Qfe+fc'("fc')exp . (37) 

• Pick up at random a set of couplings J caracterising a new function node, from the a priori distri- 
bution of couplings. 

• Compute a new u-survey, Qq{u) as 

QoM^C^o j dgdhP^{g)P2{h)5{u-U3{g,h))e^^{y[w3{g,h)-\g\-\h\]) , (38) 

where Cq is a normalisation constant insuring that Q{u) has an integral equal to one. 

This iteration defines a stochastic process in the space of u-surveys, which in turn defines a flow for 
Q{Q{u)), of which we would like to compute the fix point. Following ||l9|, this in done in practice by 
a population dynamics algorithm: one uses a representative population of M u-surveys from which the 
various Qe{u),£ € {1, ...,k + k'} used in the iteration are extracted. After Qo{u) has been computed, 
one of the u-surveys in the population, chosen randomly, is erased and substituted by Qo{u)- After 
some transient, this population dynamics algorithm generates sets of u-surveys which are sampled with a 
frequency proportional to the seeked Q{Q{u)). 

The point of this stochastic process approach is to avoid trying to write explicitely the complicated 
functional equation satisfied by Q{Q{u)). This is one crucial place where the cavity method turns out 
to be superior to the replica method: with replicas one performs the average over disorder from the 
beginning, and one is forced to work directly with the functional Q{Q{u)) |2^. As this is very difficult, 
people have thus been constrained to look for approximate solutions of Q{Q{u)) where the functional is 
taken in a simple subspace, allowing for some explicit computations to be done. 



5. Computing the energy and the complexity 

Here we show how to generalize the computation of the energy of sect. VI A~^ to the IRSB case. 
When adding one site 0, connected through k function nodes to 2k sites the energy shift in one given 
state is: 



'^-E' = X! i^^3t{hi,gi) + \hi\ + \gi\) - \^u3,{ht,gi)\ 



(39) 



where h(,g( are the incoming fields onto the function node number i. Let us call Pi{he) and P^{ge) the 
corresponding field distributions (the h-surveys). They induce a probability distribution Po{dE) of the 
energy change (^: 



[] [dhdgi Pi{hi)P',{gi)] sisE + Y. [-^^^{hi.gi) -\hi\- \gi] + | ^ uj,(/i,, 
i=i \ 1=1 1=1 



91) 



(40) 



Let us look at the corresponding change in complexity. The new system has -I- 1 variables and M + k 
function nodes, and its number of states at energy E, exp((A^ -|- 1)E([M -I- k]/[N -t- 1], [E/N + 1])) is given 
by: 



exp 



(A + 1)1] 



M + k E 



A^+ 1 ' A + 1 



= j d{5E) Poi6E)exp 



NY. 



M E-SE 
iV' N 



(41) 



This expression depends on the precise spin which has been added through the choice of the distributions 
Pe and P^', and of couplings J^, which appear in (p9|) . As one expects S(a, e) to be self-averaging, one 
must average the logarithm of the expressions in (|4l])over the iteration of population dynamics algorithm. 
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We denote this averaging by an overline. As SE is finite, one can expand in 5E/N in the thermodynamic 
hmit, to get (caUing as always a = M/N and e — E/N): 



E(a, e) - + 2a— = log U d{5E) Po{SE) cxp [-ySE] 



(42) 



where we used the fact that k = 3a. As in the RS case of (p5[), the derivative dY^/da can be computed 
by adding one function node to the system. For a generic function node a connected to the sites 1,2,3 
and with interaction coupling J, the probability distribution Pa{6E) of the energy change is: 



PaiSE) 



dhidh2dh3Pi{hi)P2{h2)P3{h3) 



SE- min {-hiai - h2a2 - h^a^ + ej{ai,a2,a3)} ~ {\hi\ 

(Ti,(T2,0-3 



(43) 



Let us look at the corresponding change in complexity. The new system has N variables and M + I 
fmiction nodes, so that: 



exp 



TVS 



M + 1 E 



N ' N 

After averaging over the iterations, this gives 



diSE) PaiSE) exp 



M E-SE 
iV' N 



II = log d{SE) PaiSE) exp [~ySE] 



(44) 



(45) 



Combining the two expressions ( p2| , [45| ) , it turns out that the quantity which is computed naturally in this 
scheme is the Legendre transform $(?/), with respect to the energy density e, of the complexity function 
i;(a,e) p9|, This 'zero temperature free energy' is defined precisely as: 



E(a,e) -ye = -j/$(y) 
and it can be computed from the population dynamics as: 

Hy) = $1(2/) - 2a$2(2/) 

$i(y) - 



y 



d7 



-ilog diSE) PoiSE)exp[-ySE]^ 
$2(2/) = -^logj^y" diSE) PaiSE) exp [-ySE[ 



(46) 



(47) 



(48) 



where Pq and Pa are given in 



I- 



Technically, it turns out that ^lyy) is more easily computed through the normalisation of the u-surveys. 
When we compute a u-survey Qoiu) as in (|3^), we can memorize the corresponding normalisation con- 
stant Cq. Picking up k u-surveys Qi at random in the population, and calling the corresponding 
normalizations, one gets: 



*,(.)^-iiog /n 



du 



Qeiue) 



exp 



yl'^uel 

. 1=1 



(49) 



VII. THE CAVITY METHOD APPLIED TO THE RANDOM 3-SAT PROBLEM 

A. Known results on the phase diagram 

Considering the random 3sat problem where the graph is generated at random and the various cou- 
plings take values ±1 with probability 1/2, numerical experiments have provided a detailed study of the 
probability PNia,K) that a given F including M = aN clauses be satisfiable. For large sizes, there 
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appears a remarkable behaviour : P seems to reach unity for a < ac{K) and vanishes for a > ac{K) 
|l3| . Such an abrupt threshold behaviour, separating a SAT phase from an UNSAT one, has indeed been 
rigorously confirmed for 2-SAT, which is in P, with ac(2) = 1 || ||. For larger K > 3, K-SAT is 
NP-complete and much less is known. The existence of a sharp transition has not been rigorously proven 
yet but estimates of the thresholds have been found. The present best numerical estimate for ac ed K — S 



is 4.26 and the rigorous bounds are|35, |3g, ^ ^ 3.26 < ac < 4.506 , while previous statistical 
mechanics analysis using the replica method, has found ac(3) ~ 4.48 and ac{3) '-^ 4.396 |Q in the 
framework of variational approximations. 

The interest in random K-SAT arises from the fact that it has been observed numerically that hard 
random instances are generated when the problems are critically constrained, i.e. close to the SAT/UNSAT 
phase boundary ^ . The study of such hard instances represent a theoretical challenge towards a 
concrete understanding of complexity and the analysis of algorithms |]l5|. Moreover, hard random instances 
are also a test-bed for the optimization of heuristic (incomplete) search procedures, which are widely used 
in practice. 

B. The cavity analysis with one state 

In the 3sat problem, the energy of a function node, as given by (^, is 

Ea = (1 + Jfai) (1 + J2V2) (1 + J^da) /4 

and leads to: 

W3Ah2,h3) = \h2\ + \h\-eiJ^h2)0{J^h-i) 

U3Sh2M) = -Jf0{JSh2)9{J^h3) , (50) 

where the function 9{x) is defined as 

e{x) = 1 if a; > ; e{x) =0 if a; < (51) 
Let us consider the cavity iteration scheme , with in the hypothesis of there being a single state. We thus 



use the general formalism presented in sect. |VI A| . In the case of 3sat, (^0|) shows that the cavity-bias on 
a given edge Ua^j takes values either in 0, 1 if the corresponding coupling is negative, otherwise it takes 
values in —1,0. Therefore the cavity- fields are integers. (This is the reason for using the unusual factor 2 
for each violated clause in the definition (^ of the energy). Because the function uj^ is an odd function 
of one of the couplings in Ja, and these couplings are random variables taking values ±1 with probability 
1/2, the distribution Q{u) of cavity-biases must be of the form: 

Qiu) = coSiu) + ii^£2l(5(u - 1) + Siu + 1)). (52) 

Plugging this expression into the self-consistency equations ( ^l|) leads to a relation between cq and the 
weight po = V{h = 0): 

= E /3a(fc) E ( 2q ) ^o'"''(^^)'' ('^q)= cxp(-3a(l - co))/o(3a(l - co)) , (53) 

k q=0 ^ ' ^ ^ 

where Iq is the Bessel function and [k/2] is the integer part of fc/2. Let us now compute cq. From ( ^0| ) 
we find that a cavity-bias vanishes whenever at least one of the incoming fields (/i2 or /13) is zero or has 
a sign opposite to the corresponding coupling. This shows that 

CO = 1 - Prob[(/i < 0) n (g < 0)] = 1 - {^^f (54) 



We obtain a closed set of equations ( |54|j53| ) which is easily solved. The distribution P{h) of cavity fields 
is then given by V{h) = "^rPr^i^ ~ where the weights Pr are equal to: 

Pr ^ V{h = r)= Vih = -r) = ± fMk) E ( 2/+ . ) cf'^-^i'-^r^^^ ( ' ) , (55) 

k=r g=o ^ ^ ^ ^ 
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and the energy is computed from ( p3| , p5| ). 

A non-trivial solution exists for a > 4.667, with a ground state energy that becomes positive at a = 5.18. 
The prediction of this hypothesis assuming a single state is a paramagnetic SAT phase with cq = po = ^, 
and energy Eq ^ for a < 5.18, and a frozen UNSAT glassy phase with cq < 1 and £'o > for a > 5.18. 
It is known ^, || that this solution is wrong both quantitatively (the location of the transition point) and 
qualitatively (the structure of the order parameter). 

The true transition is much more sophisticated, and the many state formalism corresponding to IRSB 
is needed to unveil its structure. 



C. The cavity analysis with many states 



We introduce as before the h-surveys and u-surveys, and we use the population dynamics algorithm 



defined in sect.VIBl, It turns out to be more convenient to work only in terms of u-surveys. The 
algorithm computes at each iteration a new u-survey Qq{u) by taking k + k' u-surveys Qi{u), Qk+k' (u) 
in the population through: 

Qo{u) ^ Cq J J duiQi{ui)...dukQk{uk)dviQk+i{vi)...dvk'Qk+k'{vk') 

S{u- uj{ui + ... + Uk, vi + ... + Vk')) exp[?/wj(Mi + ... + u^, vi + ... + Vk')] . (56) 



Here, k and k' are two iid random numbers taken from the Poisson distribution f^aik) defined in (|12[), 
and J denotes a set of three iid random numbers Ji, J2, J3, each being equal to ±1 with probability 1/2. 

The functions it and w are defined in (^0|). A u-survey always takes the simple form Qo{u) = (1 — 
c)d{u) -\- cS{u + Ji), it is thus a probability distribution which can be characterized by a single number c, 
and therefore the iteration of the population dynamics is easily done numerically. 



D. Solution of the self-consistency equations 

Apart from the RS solution with y — Q and Qj{u) = 5{u — Uj), where the Uj are iid taken from a 
distribution Q{u), the numerical solution finds one other solution in the region a >^ 4. Generically, the 
u-surveys found can be of three types: 

{S{u) ('trivial' or type a) 

il-7^,e~y)6{u)+7j,e-ySiu^l) (type 6+) (57) 

(1 - rj,e-y)6{u) + r,,e-yS{u + 1) (type &_) 

The arbitrary factor e~y has been introduced for convenience because numerical simulations show that 
the weight in it = ±1 of the non trivial u-surveys scale proportionally to e~y at large y. 

The statistical symmetry of the problem due to the fact that the couplings take values ±1 with prob- 
ability 1/2 implies that the probability of finding in the population a "&+" message is equal to that of 
finding a message. We call (1 — t)/2 these probabilities, and t the probability of finding a type a 

message. Non trivial messages are fully characterized by the distribution p{r]) of the i]i variables. For any 
?/, the full solution of the problem is given by the value of t and of the function p{ri). It can be obtained 
numerically by averaging over many iteration of the population dynamics. We shall now show how one 
can get some analytic control in the large y limit. 



1. Existence of non trivial u-surveys 



Looking at the iteration equation (pq), the only way one can obtain a trivial u-survey Qa{u) — d{u) is 



when either ^2 X]i=i — ^ ^^'^ whole integration domain, or J3 



< in the whole integration 



domain, or both. The probability to have J2 < in the whole integration domain is: 



Prob 



J2 ^ < 



i=l 



k=0 

exp 



2 



'2 + 2- 



l (l-t 



,k-2 



-3a- 



(58) 
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The first term t*^ corresponds to all Qi{u), i G {l,...,k}, being of type a. The second term corresponds 
to: (one of the Qi{u) of type 6+, all the other ones of type a, and J2 — —1) or (one of the Qi{u) of type 
6_, all the other ones of type a, and J2 ~ +1). The rest of the series is easily obtained similarly. 
The probability t of having a trivial u-survey is thus: 



t = 1 





k 




k' 


^1 - Prob 


J2 ^ < 


^ ^1-Prob 


J3Y.V,<0 











1 — exp 



-3a- 



1 - t 



(59) 



For a small the only solution is the paramagnetic one i = 1. A solution different from i = 1 appears above 
ao = 1.63694 (which corresponds to t — 0.04883). (In fact there appears a pair of nontrivial solutions, 
the relevant one is the one with lowest t, as can be seen from its behaviour at large a, and also checked 
in the population dynamics algorithm). It is interesting to observe that this threshold coincides with 
the point where the so called "unit clause literal" algorithm ceases to converge. At large values of 
a, the fraction of type a u-surveys becomes very small, e.g. t — 0.0051,0.0011,0.00025, at a = 4, 5,6 
respectively. 



2. Expansion at large y: location of the phase transition 



Here we shall show that, at large y, the scaling in e^^ of the weights at m = ±1 for nontrivial u-surveys 
is indeed consistent with the iteration of the population dynamics (|5^), and we deduce from this iteration 
a self-consistent equation for p(ri). 

Let us study one given iteration of (|5^). The probability of having k cavity-biases Qi, Q^, among 
which m are of type 6+, n oi type and k — m — n oi type a is: 



Ck,m,n — fsa^k) 



k^ 1 _ + 
^k — ni—n^_ _yn+n 

m\n\(k — m — ny. 2 



(60) 



Without loss of generality, we can assume that the m + n non-trivial u-surveys are Qi, Qm+m and the 
couplings are Ji = J2 = J3 = 1 for this iteration. We denote by rj the vector of weights rj = (771, ...^rjm+n)- 
Since ?ij(ui + ... + Ufc, ui + ... + t;/c') = 1 if and only if ( Ji Ui > 0) and ( J2 J2j '"j > 0)i the new u-survey 
Qoiu) depends only on the probabilities 



fM{r^) = Proh[ui + ... + Uk = q] 



(61) 



and is given by: 

Qo{u) = Co{S{u~l) 



q=l q' = l 



{m' ,n') , y{\q\ + \q'\-l) 







lEE+EE+E E 4"-)(,)4"' "')(^')e-(i^i+i''i) 

i <J— 1 q' — — k q—~k q' — 1 q— — k q' — — k I 

= Co{6{u-l)Aoe-y + 5{u)Bo] , 
Iterating the population dynamics, one finds that p(ri) satisfies the equation 

k—l k — ra oo k' — 1 k' —r 



(62) 



Y~t £ ^ 51 E E E 



Ck'. 



'=1 



An 



(63) 



This is an exact self consistent equation for the distribution p(r\) . It can be simplified at large y , and we 
show in the appendix A how to write a more tractable self-consistent equation in this limit. This equation 
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is best written in terms of the probabiHty distribution function S{4>) of the variable tf) = log(l + i]). One 
finds that S{(p) satisfies: 



(64) 
(65) 



S{(f>) = I dx+dx-dx'+dx'^A{x+)B{x-)A{x'+)B{x'^) 

(e^+ - l)(e^- - 1) 



log 



(e^+ - l)e^- + (e^+ - l)e^- + e^-e^- 



where A(x) and B{x) are two probability distributions related to S{(j)) through its Fourier transform 
S{q) = j d(l)ex:p{iq(j))S{(j)): 



A{x) 
B{x) 



1 



g3Q(l-t)/2 _ 
1 



|ie--iexp 



53Q(l-t)/2 / 27r 



— e ^ exp 



^(1-^)^(9) 
1 



- 1 



3a(t-l + -(l-i)^(g) 



(66) 



Equations (|65|,|66|) are simple enough to be solved numerically to high accuracy. 

Once S{(j)) (and therefore pirfj) is known, one can deduce the value of the zero temperature free energy 
$(y). To leading order at large y, we show in the appendix A that <&(?/) = ^/j/, with: 



dxdzB{x)B{z) log (e^ + 
-3a J Y\dxtdytB{xt)B{zi)log 

i—l 
3 

+2a / Y\_d-XidyiB{xi)Bizi)log 



1) 

2 2 

J|(e-.+e^--l)-n(e^'-l) 

3 3 

n(e--+e-^-l)-n(e---l) 



(67) 



In order to compute the S{(l)),A{x), B{x) functions, solutions of eqs ( |65| , |66|) , we have used the fact that 
these functions are probability distributions, and we have developed a population dynamics algorithm 
which follows a population of N variables (f)i, ...,(/) at, for a given value of a. 

• 1) Compute t, the solution of t = 1 — (1 — exp[3Q!(t — l)/2])^. 

• 2) Initialize the (f>j as iid random positive variables, for instance with an exponential distribution of 
width 1. Initialize the 'time variable' r = 1. 

• 3) Upgrade the time t r + 1. 

• 4) Generate an integer A; > 1 with the distribution j'^/kl l/(exp(7) — 1), where 7 = 3a(l — i)/2. Pick 
up k integers ii, .., ik at random in {1, N}, and compute the sum (/>ij + ... + 4>i^, . The distribution 
of the variable x+{t) ^ (pi^ + ... + ipi^. is A(x_|_), related to S{(j)) through (|66[). 

• 5) Generate a random variable x_(t), which will be distributed according to i?(x_), as follows: with 
probability exp(— 7), one takes x^{t) — 0; with probability 1 — exp(— 7), one repeats the procedure 
of 4) and calls the output x^{t) — (pi^ + ... + 0^^.. 

• 6) Repeat the steps 4) and 5) to generate two other variables x'j^{t) and x'_{t). 



7) Compute x(t) = e^-(^V(e''+^^^ " 1): x'{t) 



/(e-+(^) - 1), and <^o(t) - log[l + 1/(x(t) + 



X'(r)+X(r)x'(r))]. 

• 7) Replace one randomly chosen variable in the population, by the new value (f>oiT). 

The steps 3) to 7) must be repeated a large number of times, say T. One can compute the average of any 
function of 6 as: 



r 4 ^ 

j d(l>S{4>)f{cf) = — f(Mr)) 



(68) 



T=r/4+i 
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FIG. 3: Free energy ^(j/) versus the reweighting parameter y for the random 3sat problem at a = 4.2 

(One forgets the first T/4 iterations in order to reach a stationnary regime). The integrals involving a 
variable x distributed according to B{x) can be estimated by summing functions of x-. For instance, the 
first term in the zero temperature free energy ( |67| ) is evaluated as 

- dxdzB{x)B (z) log (e^' + e'-l) = — ^ log (e'-^'^ + e"-^")' - l) , (69) 

T=T/4+l 

and the two other terms are computed similarly, using the x_, x'_ values from two (resp. three) successive 
time steps. 

In practice, we used values of ^ 10000 and T ~ lOOOA^ which was enough to reach the precision 
given below on the results. 

For a < ad = 3.921, the algorithm converges towards the solution 0i — ... = 4>n = 0. This is the 
paramagnetic solution where all the u-surveys are trivial. 

For a > ad — 3.921, we find a new solution with a non-trivial distribution S^cj)). Computing the leading 
large y behaviour of the zero temperature free energy function, $(y) ^ ^'/y, on this solution using (|67|), 
we find that ^' is negative for a < ac = 4.267, while it is positive for a > ac (in the first version of this 
paper, we had reported the shghtly too small value ac — 4.256 quoted in||l^, because of a poor random 
number generator used in the population dynamics). 



E. Phase diagram of the random 3sat problem 

The previous large y analysis is in agreement with the direct numerical iteration of the population 



dynamics of sect. VI B 4, but it allows to get a much more precise determination of the thresholds. The 
results of this section have been obtained through the combined use of the numerical and analytical 
method. 

For a < ad, the system is in the SAT phase, the solution is paramagnetic, it is easy to find a solution. 
Note that, although the t equation (^^ in principle allows for the existence of nontrivial knowledges above 
a ~ 1.64, we have not found such a solution and the only one which remains is the paramagnetic one with 
t = 1. 

For ad < a < a^ we find a monotonously increasing ^(y) function, which reaches its maximum $ ^ 
at y ^ cxD (see fig. ||). 
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FIG. 4: Complexity E versus the fraction of violated clauses e = e/2 for the random 3sat problem at a = 4.2, 
obtained from the Legendre transform of fig.^ 



Fitting ^{y) by a function X^fci 6xp(— y^)/?; with p £ {1, 2, 3} gives a good and stable fit, from which 
the Legendre transform ( |4^ ) is easily done. This allows to reconstruct the complexity curve S(a, e) (see 
figj^). It is found to be finite down to e = 0. This implies that there is an exponentially large (in A'') 
number of states with zero energy density. 

The complexity of these ground states is plotted in fig||. We call this phase the Hard-SAT-Phase (HSP), 
since in this regime the typical sample is SAT, but the proliferation of states (most of which have strictly 
positive energy) makes it difficult to find a solution. 

For Qfc < a, the function <&(?/) has a maximum at a finite value y* , and $(y*) > 0. The complexity 
curve S(a, e) starts at a positive energy density eq = ^{u*) (see fig.|l). This energy density eo is the 
minimal number of violated clauses per variable which will be found in almost all samples at large N. It 
is plotted in figj|. We are in the UNSAT phase. 

This figure also shows the value of e where the complexity curve S(q;, e) is maximum gives the energy 
where there exists the largest number of states. This is the 'threshold' energy density eth where a simple 
zero temperature Metropolis algorithm (ZTMA) will be trapped. This implies that ZTMA should find 
satisfying assignments only for a < ao, in agreement with the numerical results of These predictions 
can be tested most clearly through their generalization to single instances which we discuss in the next 
section. 

Previous statistical mechanics attempts at finding this phase diagram culminated in powerful variational 
approximations using the replica method, see ref. ]4C| ] for the first results and ref. J^], which predicted 
approximate values for the SAT/UNSAT threshold - ac 4.48 in the case of ref. [|40| and ac ^ 4.39 in 
p2[ " with an intermediate phase appearing above as — 3.96 | 

The difference between the variational results and our cavity result is both quantitative and qualitative: 
in ref. the predicted nature of the intermediate phase is different with respect to ours while in ref. 
1^ the structure of the order parameter is oversimplified. In the present approach (as well as in |^) we 
work directly at zero temperature [T — 0), which has the advantage that we do not need to study the 
subtle question of the limit T 0. The reason why this limit is subtle is due to the fact that some of 
the local fields, at low temperatures, vanish linearly in T, and thus contribute to the local magnetization 
m = tanh(/377); we call these fields evanescent fields. The local magnetization at T = is zero for a zero 
field, it is equal to 1 for a finite field, and it takes an intermediate value m e] — 1, 1[ for an evanescent field. 
The variational approach of ]4C| ] focuses onto evanescent fields, and finds a continuous phase transition 
at as — 3.96 where the evanescent fields in different states start to cluster. However as these are all 
evanescent fields, this means that the corresponding local magnetizations, in a given state, are not frozen 
to ±1 but take some intermediate value, even in the T — + limit. In our T — cavity approach (as well as 



40 and Ud ^ 3.94 with the Ansatz of Ga 
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FIG. 5: The phase diagram of the random 3sat problem. Plotted is eg = eo/2 (full line), the number of violated 
clauses per variable, versus the control parameter a which is the number of clauses per variable. The sat-unsat 
transition occurs at a — ^ 4.267. The dashed line is e[i^ — tthl'i, the threshold energy (divided by two) per 
variable, where local algorithms get trapped. The dotted line is the complexity E of satisfiable states, equal to 
times the natural logarithm of their number. 



in iji^ ), the HSP corresponds to a discontinuous transition at zero temperature, involving fields which are 
not evanescent, but are of order one [ p] . This means that, in a given state, a finite number of local fields 
are non-zero integers, giving rise to magnetizations ±1, as one could expect at zero temperature. This 
phenomenology cannot be found by considering evanescent fields. Its study with replicas would require 
using a more complicated Ansatz. 

Note that our approach working directly at T = also has its limitations, for instance we are unable 
to determine precisely the self overlap (or the typical radius) of a state, or its internal entropy, precisely 
because we do not control the evanescent fields. We leave for future work the delicate study of the small 
T region of the phase diagram. Let us just notice here that a population dynamics study of this region 
with the IRSB finite temperature Ansatz of shows that the distributions of local fields tend to peak 
on integers when the temperature goes to zero in the HSP . This is a strong argument in favor of the 
exactness of this IRSB solution (i.e. the fact that we do not need to go to higher order RSB), as argued 




VIII. SURVEY PROPAGATION: CONFIGURATION SPACE ANALYSIS ON A SINGLE 

INSTANCE 



The analysis of the iterative equations for the probability distributions of messages of the previous 
sections was done by using a population dynamics algorithm, which performs an average over the under- 
lying random factor-graphs. At each step of the iteration, a random choice of coupling constants as well 
as neighborings nodes is performed with the proper probability distribution. While such an averaging 
step is central if one wants to estimate typical properties, the iterative equations make perfect sense on 
a single specific instance. The order parameters arising from the cavity equations, namely the u-surveys, 
are histograms of probability distributions of cavity biases. They determine the bias of each spin in all 
metastable states of a given energy density for a given instance of the underlying factor-graph. This is a 
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very important piece of information which can be exploited to study specific problems and to invent new 
algorithms. 

In the large N limit, we expect that the cavity assumptions hold for locally tree-like factor graphs and 
we may use u-survey propagation to have access to the properties of optimal states of minimum energy. 
Here we shall develop one such application for random 3-sat which is a concrete world-wide benchmark 
for search algorithms. However the idea of exploiting the information on optimal states carried by the 
functional order parameter is rather general and we expect algorithmic applications in different fields. 

Whenever the factor graph representing the problem does not lead to clustering within states, in practice 
whenever loops are short enough, one should think to the present approach as a first step of a sequence 
of possible approximations. 

As is well known, the so called Cluster Variation Method provides a systematic scheme that can 
be adopted to improve the approximate results given by the cavity approach W6l. The latter corresponds 
to the so called Bethe approximation which is the first step in the cluster variation scheme. Our present 
approach deals the Bethe approximation in a frustrated case. While a great deal of work has been done 
concerning higher order cluster approximations for simple models, the corresponding analysis for frustrated 
systems such as spin-glasses or hard-combinatorial problems over non locally tree-like graphs is largely 
unexplored. 



A. The Survey Propagation (SP) algorithm 



In the ordinary sum- product algorithms ||21| as described in sect. HI, the messages arriving at a node 
are added up and then sent to a function node. Next, the function node transforms all input signals into 
a new message which is sent to the descendant variable node. At each time step, on the links of the factor 
graph there are signals traveling, just like in a communication network. SP works with the same principle, 
but now the messages traveling along the links of the factor graph are u-surveys of usual messages over the 
various possible states of the system at a given value of the energy (or rather, in practice, at a given value 
of the reweighting parameter y). Of course this higher level of description is useful only when there are 
many states, which will be typically the case in hard optimization problems. One a-priori drawback of the 
approach is that the messages are complicated, being functional probability distributions: a cavity-bias 
is already a parametrization of a probability distribution (which turns out to be parametrized by a single 
variable in our case of binary spins, but could be more complicated in general); Here the u-surveys are 
probability distribution functions of these cavity-biases. In cases like Ksat at T = in which the standard 
messages can take only few values (say r), a u-survey is given by the probabilities of these values, i.e. by 
r — 1 real numbers and the SP can be implemented easily. This is one big advantage of working directly 
at r = 0, but we believe that the SP method could also be used more generally at finite temperature or 
with continuous variables, by using well adapted parametrizations of the cavity-biases. 

SP is defined for one given value of the reweighting parameter y and one given instance, with N variable 
nodes and M function nodes. Its basic ingredients are the u-surveys. Each edge a ^ j from a function 
node to a variable node j carries a u-survey Qa^j(u). The algorithm finds these u-surveys by a message 
passing procedure detailed below, and finds simultaneously all the h-surveys Pi—,a{h). Once these are 
known, one can compute the so called local field distributions and the 'zero temperature free energy' for 
this instance . The local field distribution Pi{H) on a variable node i is the distribution, over all states 
selected by the reweigthing parameter y, of the total local field H acting on spin ct^ (see (|[)). It is given 
by 



P.{H) = CJ n dUaQa^^{Ua)5 [h - ^ uJexpL ^ Ua\\ 
a£V(i) \ a£V(i) ) \ aGV(i) / 



(70) 



Ci being the normalization constant. 

We show in appendix B that the zero temperature free energy $(y) density of this sample can be 
computed as a sum of contributions ^[{y) for each function node a, corrected by the contributions ^iiy) 
from each variable node i, weighted by a factor — 1, where rii is the connectivity of variable node i: 



/M N \ 

*(y) = ]vU^*^'(2/)-E^ny)K-i) , 

\a=l i=l J 



(71) 
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where: 




beV(i)-a 



exp 



-y mm 

{ai,ieVia)} 



beV(i)-a 



--log< / Y\_ dUaQa~,i{Ua)exp iy\ E "all 1'=^- 
y [-^ aeV{i) \ aeVir) J J ^ 



= — log(a) 



(72) 



The form (JTlj) is the famiUar one for the free energy in the Bethe approximation |24j , and indeed one gets 
back the usual result using (fz^ ) in the y ^ limit. The generalization to y ^ given in d?!) adds the 
effect of the reweighting terms due to level crossings. The origin of these terms is exactly the same as in 



sect. VI 



Let us now explain how SP works. We start with a general presentation of the algorithm, which applies 
to any optimization problem involving binary variables, characterized by a given factor graph. Some 
details of the implementation for the 3sat problem will be given below. 

(0) All the u-survcys Qa^iiu) are initialized randomly. 

(1) Function nodes are selected sequentially at random; for each such node a, we update the u-surveys 
as follows (see fig. 



(1.1) for each variable node i connected to the selected function node a, we compute the h-survey 
Pi^a{h) as a reweighted convolution, see fig. (^, 

Pr^a{h) ^C^-^a j dui- ..dukQbi->i{ui) .. .Qb^^i{uk)5 - 'Y^u^ exp ^ylE""!^ C^^) 

(1.2) successively, the u-surveys on all edges a i connected to a are updated using these h-surveys: 



Qa^i{u) = Ca^i J dgdhPj^a{g)Pe^aih)S{u~u,j{g,h))exp{y[wj{g,h)-\g\-\h\]) (74) 

{Ci^a , Ca^i are normalization constants) . 

(2) The iterative process of step (1) continues until convergence is reached. If the process converges, 
the corresponding stable set of u-surveys are used to compute the N local field distributions using 
([70|), and the zero temperature free energy <&(?/) given in ( [7l| , [7^ ). 

The above procedure can be repeated for different values of the reweighting y so that the complexity 
S(2;) = d^{y)/d{l/y) and the energy density e{y) = d{y^{y))/dy of states can be estimated. The 
parametric plot of S(?/) versus e{y), varying y, gives the complexity S(e) of states of energy E = Ne. 

When it converges, SP allows to get an order parameter (the set of all the surveys), a zero temperature 
free energy density $(y), and a complexity curve S(e) for one given instance. What is the meaning of 
these quantities in general is an open question. A given instance has a finite value of N , and therefore 
the notion of 'state' is not easy to define. Roughly speaking one can think that for large N, there might 
exist some effective 'finite N states', such that the number of spins to flip in order to reach one state 
from another one is large, leading to a separation of scales of the number of spins involved between the 
intra-state moves and the inter-states moves. Such a situation would generally be difficult to handle for 
search algorithms, and this is where SP could be quite useful. In order to get a first understanding of 
these questions, we have experimented SP on single instances of the random 3sat problem for large values 
of iV. 
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FIG. 6: Function node a and its neighboring graph 



B. The case of random 3sat 



We consider one instance of the 3sat problem, chosen randomly as in sect. with energy 

M 



n 

Q=ijey(a) 



1 + Ja^jO'j 



(75) 



For 3sat the cavity-biases on a given link a ^ j takes values Ua^j = 0, —Ja-ij', The corresponding survey 
is Qa^j{u) — Ca^j5{u) + (1 — Ca-tj)S{u — Ja^j)- The full set of u-surveys is characterized by the 3M 
numbers Ca^j which are updated according to the SP algorithm described above, until convergence. The 
results of our numerical experiments are the following. 



1. The paramagnetic phase 



For a < ad, SP converges toward the trivial paramagnetic solution Qa^iiu) = S{u), for all a — > i edges. 
Local field distributions are also trivial, Pi{H) — 6(H) Vi, and no information can be gained on the fine 
structure of ground states. In this region, there is a single state and the statistical properties of the zero 
energy configurations are totally driven by its entropy. A different formulation of the cavity approach, in 
which the proper /? — > oo limit is taken and the evanescent fields are computed, could reveal some finer 
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FIG. 7: Extensive zero temperature free energy N^{y) versus reweighting y for two specific instances of size 
N = 10000 with M = 42000 and M = 45000 clauses respectively, (a = 4.2 and a = 4.5). 

information for this paramagnetic phase, which however is known to be trivial from the algorithmic point 
of view. 

2. The intermediate phase 

For ac > a > a^, that is in the glassy region, the random sequential updating of the iterative process 
converges to a unique non-trivial solution, provided y is large enough. In practice, we start from y large, 
like e.g. y = 6 (remember that the corrections to the y oo limit are exponentially small), run SP, and 
after finding a solution for the u-surveys we decrease y (e.g. y y ~ .2) and rerun SP using the previous 
u-surveys as a starting configuration for this new y. This speeds up the convergence. Below some value of 
y the non-trivial solution disappears abruptly and the algorithm converges to the paramagnetic solution. 

In this region of a, the solution space as well as the configurations of higher energy become divided 
into an exponential number of states. To compute the complexity, we measure the free energy $(?/), see 
fig. (^, and we perform the Legendre transform numerically. 

The curve of the total complexity NY. versus the total energy Ne for one sample of random 3-sat with 
N = 10000 and M = 42000 is given in fig. (|). One finds NY{E = 0) 34, meaning that the zero energy 
(SAT) states are predicted to be exponentially numerous, e"^^ at the leading exponential order (remember 
that each such state itself contains a large number of spin configurations |j]). The threshold states have 
an energy of approximatively 44 violated clauses and their number is predicted to be about e^^^. A cross 
check of such predictions is given by the behaviour of ZTMA which cannot cross energetic barriers. It 
can be shown that for random 3-sat zero energy moves allow to explore configurations within each state 
and therefore we expect such algorithms to get trapped in the most numerous ones (the threshold states). 
Indeed, extensive numerical simulations of ZTMA on many samples of different size (ranging from few 
hundreds to 10^) and for different values of a confirm such scenario. As a representative example, we 
report that for the sample whose compexity is plotted in fig. (^), repeated runs of ZTMA get stuck at an 
energy sharply peaked around 48 violated clauses, with a small residual dependence of the final energy on 
the simulation time (the final energy found by ZTMA shows a power law behavior on the total number 
of spin flips). 

We have checked that the functional order parameter given by the u-surveys and the local field distri- 
butions carry precise information concerning the space of solutions for one given sample. Working with 
N not too large, some SAT configurations can be found efficiently by good algorithms like e.g. walksat-35 
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FIG. 8: Extensive complexity, A'^E, versus the total number of violated clauses NE' (= NE/2), for the specific 
instance of size N = 10000, M — 42000 studied in figj^ The complexity is obtained as the Legendre transform of 
the zero temperature free energy. 

p8| , |49|| . We have collected a large set (1000) of uncorrelated SAT configurations by running this algorithm 
many times with random initial conditions. In each such configuration w, the spin Ci takes a value af , 
and we have computed, for each given site i, the average of the whole set SAT configurations oj. Next we 
have compared the above results with the predictions of SP as follows. 

We first have selected the states with minimal energy, by picking the value of y which maximizes F{y). 
Here this is y = cxd and for practical computation it was enough to choose y sufficiently large (corrections 
are exponentially small). According to eq.([70|), the field Hi — J2a "^a-ti in each state is an integer valued 
variable which can be computed from the u-surveys. The total weights 

/•oo />0^ 

w+ = / dHP,{H) ; w- ^ dHP,{H) (76) 

of Pi{H) on positive (resp. negative) integers give the fractions of zero energy states where ai is fixed to 1 
(resp. to —1). As displayed in fig. we find a remarkable agreement between the local magnetizations 
wf — w~ predicted by SP and the local magnetizations measured by averaging over the ground states found 
by the walksat algorithm. In the figure we report data for N = 10000 and M — 42000: we have devided 
the local magnetization in 30 intervals and labeled spins according to the prediction of SP. Next on such a 
partitioning of spins we have taken the average over the configurations found by walksat. (The remarkable 
agreement of numerical and SP results indirectly shows that walksat is a good uniform sampler.) 

The weight in = of Pi{H), w^=l~ wf — w" , measures the tendency of a variable to be under 
constrained: for instance, variables which belong to very few clauses have = 1. 

3. The UN SAT phase 

For a > ttc, SP predicts a positive ground state energy with zero complexity, whereas excited states 
remain exponentially numerous. Proper tuning of the reweighting, that is choosing y so that the complexity 
vanishes, allows to predict the ground state energy and to evaluate the probability distribution of effective 
fields for each variable. In this regime SP is found to converge only when the reweighting parameter is 
well chosen. For small values of y, SP converges to the paramagnetic solution or to the RS solution. 
For intermediate values of y, SP converges to the non-trivial solution whereas for larger values of y, SP 
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FIG. 9: The bias of the variables predicted by SP (with y = 8) compared with the one measured analyzing SAT 
configurations from the same sample (A^ = 10000, M = 42000) 



stops converging. The range of y values for which SP converges to the non-trivial solution is sufficient to 
determine the free energy. An example is given in fig.^. 

For large values of a, we expect multiple nested clustering phenomena to appear, that is continuous 
replica symmetry breaking[^. This scenario could be analyzed by a further generalization of SP which 
is beyond the scope of this work. 



IX. SURVEY PROPAGATION AS A SOURCE OF NEW ALGORITHMS FOR HARD 

OPTIMIZATION PROBLEMS 



The previous section has shown how SP can give rather precise answers on the structure of the space 
of configurations and the ground state energy of the random 3sat problem. Here we shall stay within this 
problem and ask the following natural question: 

Given a random 3-sat formula of size N, how can we take advantage of SP in order 
to find optimal configurations? 

If SP could predict with very high accuracy the value of the GS energy of a given formula, it could also 
predict its satisfiability. Then one could proceed in finding a satisfying assignment just by converting the 
decision algorithm into a search algorithm as follows A variable is selected and fixed to one value. We 
then use SP to evaluate the GS energy of the subproblem of size — 1 and decide whether it is still SAT 
or not. If the subproblem is SAT then we keep the assignment, otherwise the opposite value of the binary 
variable is chosen. The process is repeated untill all the variables have been exhausted (in at most 2N 
steps). If along this reduction process the subsystem becomes a paramagnet, then SP becomes ineffective 
and another search algorithm must be used on the subsystem (but for a paramagnet it is very easy to find 
the ground state). 

The above scheme, however, suffers from finite size effects and from the imprecision in the determination 
of the ground state energy, a fact which is particularly important close to ac- Moreover it does not take 
advantage of the information provided by the u-surveys. 
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FIG. 10: Effect of fixing a single balanced spin on the complete complexity curve A^E versus the number of violated 
clauses ( i.e. one half of the extensive energy) of an instance of 3sat with = 1000 and M = 4200. The difference 
in the two curves is very close to In 2. 

A. Categories of variables in one specific instances 

A somewhat coarse grained information contained in the u-surveys, once SP has reached convergence, 
is the one given by the total weights wf^ of the local field distribution which gives the fraction of states 
where the spin ai is positive (resp. negative). Having computed these weights, we may distinguish three 
reference types of spins (of course all the intermediate cases will also be present): the paramagnetic ones 
with ~ 1, the biased ones with wt ^ 1 or w~ ~ 1 and the balanced ones with wt — w~ and u;? small. 

In order to characterize the differences between these various types of variables, we have performed a 
few numerical experiments and analyzed the effect of fixing one such spin on the structure of states of 
the subproblem of size N — 1. As displayed in the complexity curve of fig.(^0|), the three types of spins 
produce different effects, consistently with the interpretation of the order parameter. Fixing a biased spin 
does not alter the structure of the states and the complexity changes smoothly. Fixing a paramagnetic 
spin has an effect only on the internal entropy of the states (which we cannot measure) but leaves the 
energy unaltered. Interestingly enough, balanced spin have an enormous effect: the most balanced ones 
produce a decrease very close to In 2 in the complexity, indeed half of the states are eliminated by fixing 
one single balanced variable!. 

B. Survey inspired decimation (SID) algorithm 

One strategy for using this information in order to produce an optimization algorithm is to fix as many 
variables as possible without altering the ground state energy, evaluated step by step as the size of the 
problem decreases. Eventually, either all variables have been fixed or (more likely) the remaining variables 
turn out to be paramagnetic (i.e. Pi{H) = S{H), Vi), in which case a simple search process can be run to 
find the complete ground state configuration. 

A straightforward implementation of the above ideas provides a simple algorithm that can be used to 
find solutions to random 3sat in the hard region a G [ad, ad- We do not expect this implementation to 
be the most efficient one, in that no particular strategy has been worked out to optimize the decimation 
process. The scope of this first implementation consists in showing the potentiality of the novel algorithmic 
scheme and we leave for future work the design of optimized versions of the algorithm or applications in 
different contexts. 
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The overall idea underlying the search process is rather simple. At each time step a single variable is 
fixed according to the outcome of SP and the effect of such fixing is used to simplify the problem. The size 
of the problem reduces from Nt to Nf—l — St, where St is the number of variables which become fixed due 
to the simplification of the problem: satisfied clauses are eliminated, unsatisfied K clauses are transformed 
into (K-1) clauses. K — 1 clauses need to be satisfied and therefore their variables are properly fixed (unit 
clause propagation) leading to further spin elimination. 

At the beginning of the process, randomly chosen balanced spins can be fixed in order to reduce the 
number of states. At each step one may compute the free energy to detect the onset of violated clauses. 
One may also evaluate the function <&(?/) to have an estimate of the complexity. Successively, biased spins 
are fixed. Whenever a paramagnetic state is found, or at some intermediate steps, a rapid search process 
like simulated annealing at a fixed cooling rate or walksat is run on the sub-system. We may end up 
either by having found a solution or by having still few violated clauses. In the latter case we may simply 
restart. The sketch of the SID algorithm is: 

0. Random initial condition for the cavity-biases 

1. Run SP and evaluate {Pi{H)} , or {wf,w~,w^}, and . 

2. Check for a paramagnetic state and in case (or at some intermediate step) run a 
fast local search process (e.g. simulated annealing or walksat). If a solution 
is found output ''SAT'' and stop. 

3. Select and fix the most biased variable (the one with the largest jw^ — w^\) and 
simplify the problem. 

4. If the problem is solved completely by unit clause propagation, then output 
''SAT'' and stop. If no contradiction is found then continue the decimation 
process on the smaller problem (go to 1.) else (if a contradiction is reached) 
restart (go to 0.) 

Extensive numerical experiments on random 3sat instances at a = 4.2 with size up to = 10^ have 
shown a remarkable efficiency of SID. While the process of fixing a single variable takes some time {0{N) 
operations) the number of assignments explored is very small. At a = 4.2 typically a single run of SID 
(i.e. with no restarts) leads to a solution. Closer to the critical a, few restarts might be necessary in order 
to find a configuration of strictly zero energy. However, at each run the typical energy found by SID is 
very close to zero, well below the energy at which simulated annealing gets stuck. A detailed description 
of the numerical experiments will be given in a furthcoming paper pT| . We just mention that the largest 
public benchmarks of random 3sat ||49| have been solved efficiently by SID. 

In fig. (|ll|) we show the evolution of the complexity under SID. For a sample of size N — 10000 at 
a — 4.2 we evaluate the complexity curve every 200 decimation steps until a paramegnetic state is reached. 
SID acts by eliminating clusters of solutions and hence reducing the complexity of the ground state down 
to the point where very few clusters remain. 



X. CONCLUSION 



We have derived here two main results. The first one concerns the phase diagram of the random Ksat 
problem, and establishes the existence of an intermediate phase where the problem is SAT but the solution 
is difficult to find because of the existence of many states. We would like to point out that the cavity 
method which we have used here is not rigorous: it relies on some hypotheses which can be true only 
for large systems and are thus difficult to prove (although this can be done in some cases \^t\). However 
the experience gained from similar problems, together with numerical results of this and previous papers, 
indicates that this solution is likely to be the correct one. If higher order replica symmetry breaking 
effects would show up, one can believe that in any case their quantitative infiuence on the results should 
be rather small. It should also be noted that the very same cavity strategy which we have used here has 
been tested on a variant of the Ksat, the XORsat problem, which can be also solved by exact methods 
p2| , |53| . In this case, the existence of the intermediate phase has been confirmed and all the predictions 
(qualitative and quantitative) from the cavity method have been checked rigorously p^ . 

The second result gives a new class of message passing algorithms for solving optimization problems 
in the regime where the proliferation of metastable states slows down a lot all the local algorithms. One 



30 




such algorithm turns out to be quite powerful on the case of random 3sat problems. Clearly a lot of 
work needs to be done in order to develop such algorithms in various contexts and test them against 
traditional strategies. A more direct derivation and understanding of the algorithm, and in particular of 
the nontrivial reweighting term, would also be welcome. 

Finally, we expect that the single sample SP method at finite temperature can become a useful tool in 
analyzing the fine structure of order parameters in disordered systems and other complex systems. 



Acknowledgements 

We thank G.Parisi, B. Selman and J.S.Yedidia for useful discussions and encouragement, and S. Mertens 
for pointing out the inaccuracy in our first estimate of ac- RZ thanks the LPTMS, Orsay, for the kind 
hospitality. 



APPENDIX A: LARGE y EXPANSION FOR THE 3-SAT PROBLEM 

We give here some details on the solution of the population dynamics equation ( |56| ) of the 3sat problem 
at large y. We start from the self-consistent equation (|6^) for the distribution of the rescaled weights of 
the u-surveys in u = 1, and we recall that fq"^'^\ri) is defined as the probability that ui + ... + Uk — q, 
given that Qi{ui), ...,Qm{um) are of type 6+, and Qm+iium+i)-, ■■■,Qm+n{um.+n) are of type 6_. The 
quantities Aq, Bq are expressed in terms of three numbers g^^-, g^, go'. 

k -1 

/('»•") (,7)e''l'l ; g_(r;) ^ J] /(™^")(,7)e^l^l ; gM ^ flr''\v) (Al) 

g=l q=—k 

as: 



^0 = 9+{v)9+{v') 

Bo = (?+(??)[g_(7?')+^?o(?7')]+3+(^')[5-(^?)+5o(^/)] + [5-('/)+^?o(^y')][ff-W+5o(??)] (A2) 
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The next step consists in introducing the joint probabiHty distribution P^"^'"\g+,g-,go) 



9=1 



(A3) 



q— — k 



Eq. (H) reads 



00 k k—m 



00 k' k' — m' 



p(wo) = E E E "^fc,™^" E E E '-'k'.ni'y / dg+dg_dgodg'^dg'_dg'Q 

fe=Om=ln=0 fc'=Om' = l ri'=0 



A + B 



(A4) 



where the coefficients Ck,m,n are given in (|60|), and A, i? are given in (p^). 

The scahng of the weights of the u-surveys in u = 1 wq with e^^ at large y is consistent with the 
self-consistency equation (|63[). In this limit we find g+,g-,go ^ and these quantities simplify to 

(up to corrections of order e~^): 



9+iv) = (l + ??i)...(l + ?7m)-l 

9-{v) = (1 +??m+l)---(l +'/m+«) - 1 

50 (?7) = 1 



(A5) 



Therefore, the three variables g-^-,g-,gQ become uncorrelated random variables in the large y limit, with 
a distribution: 



p(™'")(5+,g_,5o) - n^"'Hg+)n^-Hg-)Sigo - 1) 



where 

^^^Hg) = / d'ni...dr]ip{r]i)...p{r]i)S | g 
The equation for p(?7o) reads: 



U=l 



(A6) 



(A7) 



00 k k—m 



oc k' k' — m' 



PiVo) = YTT^E E E E E E '^k',m',n' / dg+dg-dg\.dg'_ x 

fc=Om=ln=0 fc'=Om' = l n'=0 

X f7(™)(.g+)f7(")(.g_)17("')((7V)17("')(.g^) x 



X 5[riQ 



9+9^ 



5+(l + 5^) + 5V(1 + 5-) + (1 + 5-)(l + .9-) 



(A8) 



Clearly, the function fl^^^ (g) can be seen as the £-th convolution of a certain function after an appropriate 
change of variables. One is lead to introduce the variables 4>i and x defined by: 



l0g(l + 77i) , x = log(l + 5r) 



and we call ^(^i) , T^"^\x) their probability distributions. Equation(A7) shows that 



t(")(x) = J d(t>i...d^^S{(f>i)...Si^^)6 (^x-f^^^j 



(A9) 



(AlO) 



where (f>i £ [0,ln2] and x E [0, to In 2]. 
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In order to simplify the self-consistency equation, we introduce the joint probability distribution 

k k — m 



(All) 



k m— 1 n— 



which is normalized: ^ dx+dx-R{x+,X-) = 1. 

We may show that R factorizes by introducing its Fourier Transform: Using the coefficients (|60|), the 
triple series can be resummed and expressed in terms of the Fourier transform S{q) of S{x): 



J dq^dq-R{x+,x )e 



-3a 



E 



(3a)^ 



1 



1 



exp 



exp 



3a - 1 + 
3a ( ^ - 1 



2 

1 - t 



2 

1 -t 



S{q-i 



2 

1 - t 



t + Kr^S{q^) 



Siq-) 



-S{q-) 



Rearranging the above expression and taking the inverse transformation we find for R 

R{x+,x-) = A{x+)B{x^) , 

where 



(A12) 



(A13) 



A{x+) 
B{x-) 



1 



dq. 

g3a(l-t)/2 _ 1 / "27r 



1 



dq^ 

g3a(l-t)/2 / 



g-ig_a;_g^(l-t)S(<2+) 



(A14) 



We may now write the self-consistency equation in a tractable form. Defining the variable (j>o associated 
with r/o as: 



(e^+ - l)(e^+ - 1) 



(e=^+ - l)e^- + (e^+ - l)e^- + e^-e"- 
we transform the equation for /5(77o) into an equation for S((\)q) 

S{4>o) — J dx+dx-dx'j^dx'_R{x+, X-)R{x'_^, x'_) 

(e^+ - l)(e"=+ - 1) 



(A15) 



(5 00 - log 



(e^+ - l)e^- + (e''+ - l)e^- e^-e^ 



(A16) 



This is the equation that we have used in order to solve the problem numerically. 

Let us mention however that a series of further simplifications may also be written. It is easy to verify 
that 



S{(j)o) = J dzdz'C{z)C{z')d Uo 



where 



log 



C{z) = / dx+A{x+)dx-B{x-)6 z - 



1 



- 1^ 

Moreover, if we perform in the above the change of variable 1 + z ~ e'' , defining the distribution 

DiO = Ciz)-^ , 
dC 



(A17) 



(A18) 



(A19) 
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we find: 



S{M = / dr,E{Tj)S 00 - lo, 



e'' - 1 



(A20) 



where E{ri) is the convolution of D with itself: 

E{rj) EE I dCdC'DiODiOSiC + C-V)- (A21) 

The transformations which we have written from S—fA,B~^C'—^D^E—fS can all be done using 
one dimensional integrals, Fourier transforms, and changes of variables. Hence this provides an iterative 
mapping from the function S{(f)) onto itself which can be handled efficiently numerically. Which form of 
the self-consistency equation to use is a matter of c ompu tational convenience. It turns out that, for our 
purpose, enough precision could be obtained from (A17) and we did not try to develop this alternative 
computation. 

We now proceed with the computation of the energy ^(jj) defined in (^). We shall be interested in 
evaluating it at large y, using the solution for the distribution of u-surveys that we have just found in this 
limit. 

We start with the piece It consists of two terms: 



1, 



$i(y) = — log Co logAo 

y y 



(A22) 



where the overline denotes the average over the population dynamics of sect. |VI B 4| and where Co and 
Aq are given by 



J]^ [duiQt{ui)] exp 



1=1 



^ J dgdhPi{g)P2ih)exp[ywj{g,h)] . 
The term Aq is easily written in terms of the variables g+,g^,gQ: 

Ao = J2 /i"'"^ We^'^' = 9+iv) + 9-iv) + 9o{v) 
Averaging over the population dynamics, we get for large y: 

-logAo = - V Ck,m,n / dgodg+dg^P'^"^-''\g+,go,g^)\og[go + g+ + g-] 

y y I. J 



(A23) 
(A24) 

(A25) 



- "^fe.™'" I dxdzT(^\x)T^^\z)\og{e- + -I) 



(A26) 



Treating separately the m > 1 piece which can be resummed as in (All,A13), and the m — Q piece, the 
sum over k,m,n gives 



^ Cfe,™,„r(™)(x)r(")(z) = VT^A{x)B{z) + S{x) J2 Cfc,o,„r(")(z) - six) Biz) 

k,rn^n k,n 



Putting this expression back into (A26), we finally find 

-T^gA = - / dxdzB{x)B{z) log (e^ + - 1) 

y y J 



(A27) 



(A28) 



We now turn to the second contribution, Co, to (A24); before averaging over the iteration of the 
population dynamics, we have: 



1 



k k' 

W [duiQi{ui)] [dv.jQk+]{vj)\ exp 



4=1 



k k' 

ywj ( 

i=i j=i 



(A29) 
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Without loss of generality we assume Ji = J2 = 1; we u se th e solution (|57|), and denote as before by 
m,n,m' ,n' the numbers of various u-surveys appearing in(A29): 

k — > TO type 6+, n type 6_ and k ^ m — n type a 
k' to' type &+, n' type &_ and k' — m! — ri type a 



1 



Co can then be written in terms of the g variables as: 

dgod.9+d.g-P(™-")(.g+, .go, 5-) dg'^dg\dg'_P^^' {g\,g'^, g'J) 
[g+g'^e^y + 5+ {d- + 5o] + 5^ \9- + 5o] + [ff- + 5o] [.9- + 5o]) (A30) 
For y large we can drop the term in exp(— y); averaging over the iteration of the population, we get 

dx'dz' r(™')(a;')T("')(z') log [(e^ - + (e^' - l)e^ + e^e^' 



log 



= / dxdzdx'dz' B{x)B{z)B{x')B{z') log [(e^ - l)e^' + (e^' - l)e^ + e^e^'l (A31) 



Finally we now compute the $2(2/) piece given in (p8[). For a generic clause with coupling Ji, J2,t/3, 
involving the h-surveys Pi{hi), ^2(^2), ^3(^3), we have 



exp[-y$2(y)] - j dhidh2dh^Pi{hi)P2{h2)Pz{h^)t 



-2ye{.Jtht)e{j2h2)e(.hh3) 



(A32) 



We write as before Pi{hi) = J Y\eLi[dueQi{ui)], and suppose that this set of ki u-surveys contains mi 
u-surveys of type 6+, ni of type b- and fci — toi — ni of type a, characterized by the weights rji. The 
same decomposition is done for P2(^2) (resp P3(/i3)), where the numbers of u-surveys of various types 
are 7712 , n2 , fcg — TO2 — n2 and the weights are ry2(resp. 1713,113, k^ — TO3 — 77.3,773). For this iteration, the 
expression ( A32| ) for $2(2/) can be reexpressed as: 

exph2/$2(2/)] = /ir'"'H'?i)4r""'H^2)4™-"^n'?3)exp[-2y^(gi)0(g2)0(g3)] 



91,92,93 
3 



n [5+ + 5- (^yi) + 50 - n [5+ (^i) 



(A33) 



i=l 



The average over the population gives: 

My) ^ J2 J2 J2 Ck,,m,.n, Q 



fei.mi.m fc2,m2,"2 fc3,™3,n3 



pj(e-.+e^--l)-n(e^'-l) 



- / JJ(ia;i(izii?(xi)i3(zi)log 



J|(e-. + e^' - 1) - n(e"- " 1) 



(A34) 



Grouping together the contributions ( A22| . A28 , A3 1 , A34 ) , we find the total zero temperature free energy 
density <I>(y) = $i(y) — 2a$2(2/) given in (37|). 



APPENDIX B: FREE ENERGY FOR ONE GIVEN SAMPLE 



Let us explain here how to compute the zero temperature free energy for one given sample. We start from 
the contribution of one given factor node a. We shall look at a somewhat large part of the graph containing 
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a (see fig.|2|). We call cti, 0-2, the three spins connected to it. One of these spins, ar is connected, beside 
a, to kr other function nodes which we call b^, ...,b^'- . The function node is connected, beside ar, to two 
other spins which we call and t^, and the cavity fields onto them are called and (see fig.|l2|). In 
the absence of the spins (Ti, (T2, (T3 and of all the function nodes the ground state energy of the system 
would be 

3 fcr 



^-«=-EE(l-9rl + l^rl) (Bl) 



r— 1 s— 1 

Adding the spins ai,(72, fa and all the function nodes 6*, the ground state energy becomes 

3 A,v 



Ef,n = min i (cti , (72 , (73 ) + V V ( min [Eb^ (cr^ , a^. , ) - g^cr^ - ft,,"r;] ) 

\^ r—1 s— 1 ^ ^ , 

= min \ Ea{cri,a2,(T3) -y^ary^^ujigr,K)\ 

rr-i .rr'i.fTr, I ' * * I 



1 ,£'"2 :Cr3 
3 kr 



r—1 s—1 



-J2J2('^^(9:;K)~\9r\~\K\) (B2) 
r—1 s—1 

The zero temperature free energy shift induced by the addition of all these nodes is given by: 

where the integral is over all the , M fields, each with a probability distribution given by its h-survey. 
We can use the iteration equation (^J) in order to simplify this complicated integral through the use of 
u-surveys: 



g-y*a(y) 



3 K ^ 



k,. 

nn 77^^<Qbf---K-) 



exp 


f . j 

—y mm < 




Cri,(T2.<T3 



(B4) 



r=l s=l 



We now compute the contribution from one given variable node z = 0. We use the notations of fig.^ 
calling ar the k function nodes to which it is connected (r e {1, .., k}). The function node is connected, 
beside ctoj to the spins ar,Tr, and we call gr (resp. hr) the cavity field on cr^ (resp. r^). In the absence of 
spin (To and of the function nodes connected to it, the ground state energy of the system would be: 

k 

E^mt^-Y,{\gr\ + \hr\) (B5) 

r=l 

Adding the new spin and the function nodes a^, the ground state energy becomes: 
Efin = minV" min [£'q^(cto, cr^, r^) - grcrr - KTr] I 

(To \<T,.,Tr / 

r—1 ^ ' 

Cfc \ fc k k 

-(TQ^U3A9r,hr) \ - ^ Wj ^{gr , K) = Uj ^{gr , hr)\ ~ ^ Wj ^{gr , K) (B6) 

r—1 / r—1 r=l r=l 

The zero temperature free energy shift induced by the addition of these nodes is given by: 



where the integral is over all the gr,hr fields, each with a probability distribution given by its h-survey. 
As usual, this can be simplified by the use of u-surveys derived from the iteration eq. (|73) : 



^ '" dur 



n 



Qa^^o{Ur) 




FIG. 13: The set of nodes with which one computes the free energy shift 3>o(y)- 
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A little thought shows that, when computing the total zero temperature free energy $(?y) = ^liv) ~ 
^■(n^ — !)$"(?/), one is correctly counting each node once. In particular, in the limit of y ^ 0, $(y) 
reduces to the sum of the energy of all factor nodes, as it should. The same reasoning shows that the "1/ C" 
factors in (B4) and (B8) actually cancel, so that one can forget these normalisations for the computation 
of $(y), as was done in the text in formulae (|7l|,[72|). 
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